

full_set <- readRDS("temp/geocoded_shootings_from_gva.rds") |> 
  rename(longitude = lon,
         latitude = lat,
         id2 = id) |> 
  ungroup()

full_set <- left_join(full_set,
                      select(fips_codes, state_2l = state,
                             state = state_name) |> 
                        distinct()) |> 
  select(-state) |> 
  mutate(state_2l = tolower(state_2l))

runs <- list.dirs("raw_data/vest/vest_2016/2016")[2:length(list.dirs("raw_data/vest/vest_2016/2016"))]

ps <- read_sf(dsn = runs[1], layer = substring(list.files(runs[1])[1], 1, nchar(list.files(runs[1])[1]) - 4))

crs_hold <- st_crs(ps)

find_precinct_16 <- function(s){
  library(tigris)
  library(rgdal)
  library(sf)
  library(sp)
  library(rgeos)
  library(SearchTrees)
  library(raster)
  library(data.table)
  library(tidyverse)
  
  state <- substring(list.files(s)[1], 1, 2)
  
  ps <- read_sf(dsn = s, layer = substring(list.files(s)[1], 1, nchar(list.files(s)[1]) - 4))
  colnames(ps)[1:(ncol(ps) - 1)] <- toupper(colnames(ps))[1:(ncol(ps) - 1)]
  
  ps <- st_transform(ps, crs_hold)
  
  ps <- ps |> 
    mutate(GEOID = paste0(state, row_number()),
           share_dem = G16PREDCLI / (G16PRERTRU + G16PREDCLI)) |> 
    select(share_dem, GEOID)
  
  tryCatch(
    {
      ps <- st_make_valid(ps)
    }, error = function(e) {
      Sys.sleep(1)
      cat("Error:", e$message, "\n")
    })
  
  ps <- ps[st_is_valid(ps), ]
  
  ms <- filter(full_set, year(date) %in% c(2016, 2017),
               state_2l == state)
  
  points <- st_as_sf(SpatialPointsDataFrame(cbind(ms$longitude, ms$latitude), ms))
  
  points <- st_set_crs(points, st_crs(ps))
  
  points <- st_join(points, ps, join = st_within) |> 
    select(id2, share_dem, GEOID) |> 
    st_drop_geometry() |> 
    mutate(year = 2016) |> 
    filter(!is.na(share_dem))
  
  return(points)
  
}

find_precinct_20 <- function(s){
  library(tigris)
  library(rgdal)
  library(sf)
  library(sp)
  library(rgeos)
  library(SearchTrees)
  library(raster)
  library(data.table)
  library(tidyverse)
  
  state <- substring(list.files(s)[1], 1, 2)
  print(state)
  ## pull BG shapefiles using tigris package
  ps <- read_sf(dsn = s, layer = substring(list.files(s)[1], 1, nchar(list.files(s)[1]) - 4))
  
  colnames(ps)[1:(ncol(ps) - 1)] <- toupper(colnames(ps))[1:(ncol(ps) - 1)]
  
  ps <- st_transform(ps, crs_hold)
  
  ps <- ps |> 
    mutate(GEOID = paste0(state, row_number()),
           share_dem = G20PREDBID / (G20PRERTRU + G20PREDBID)) |> 
    select(share_dem, GEOID)
  
  tryCatch(
    {
      ps <- st_make_valid(ps)
    }, error = function(e) {
      Sys.sleep(1)
      cat("Error:", e$message, "\n")
    })
  
  ps <- ps[st_is_valid(ps), ]
  
  ms <- filter(full_set, year(date) %in% c(2020, 2021),
               state_2l == state)
  
  if(nrow(ms) > 0){
    points <- st_as_sf(SpatialPointsDataFrame(cbind(ms$longitude, ms$latitude), ms))
    
    points <- st_set_crs(points, st_crs(ps))
    
    points <- st_join(points, ps, join = st_within) |> 
      select(id2, share_dem, GEOID) |> 
      st_drop_geometry() |> 
      mutate(year = 2020) |> 
      filter(!is.na(share_dem))
    
    return(points)
  }else{
    return(filter(data.frame(id2 = NA, share_dem = NA, GEOID = NA, year = NA),
                  row_number() < 1))
  }
  
  
  
}

cl <- makeCluster(8)  
registerDoParallel(cl)

clusterExport(cl, list("find_precinct_16", "full_set", "crs_hold", "find_precinct_20"))

out <- rbindlist(parLapply(cl, runs,
                           fun = find_precinct_16))

saveRDS(out, "temp/shootings_precincts_2016.rds")

runs <- list.dirs("raw_data/vest/vest_2020")
runs <- runs[2:length(runs)]
runs <- runs[runs != "raw_data/vest/vest_2020/nj_2020_vtd_estimates"]

out <- rbindlist(parLapply(cl, runs,
                           fun = find_precinct_20))

saveRDS(out, "temp/shootings_precincts_2020.rds")

############################################################

b2p <- function(state){
  
  library(tigris)
  library(rgdal)
  library(sf)
  library(sp)
  library(rgeos)
  library(SearchTrees)
  library(raster)
  library(data.table)
  library(tidyverse)
  
  #if(!(file.exists(paste0("temp/p2b_", state, ".rds")))){
    bgs <- tigris::blocks(state = state, class = "sp", year = 2019)
    
    bgs@data <- bgs@data |> 
      mutate(across(c('INTPTLON10','INTPTLAT10'), as.numeric))
    
    
    pct <- readOGR(paste0("raw_data/vest/vest_2016/2016/", tolower(state), "_2016"),
                   paste0(tolower(state), "_2016"))
    pct <- spTransform(pct, CRS("+proj=longlat +datum=NAD83 +no_defs"))
    
    pct@data$GEOID <- paste0(state, c(1:nrow(pct@data)))
    
    pings  <- SpatialPoints(bgs@data[,c('INTPTLON10','INTPTLAT10')], proj4string = pct@proj4string)
    bgs$precinct <- over(pings, pct)$GEOID
    
    h16 <- bgs@data |> 
      mutate(state = state,
             year = 2016)
    
    if(state == "KY"){
      
      pct <- readOGR(paste0("raw_data/vest/vest_2020/", tolower(state), "_2020_vtd_estimates"),
                     paste0(tolower(state), "_2020_vtd_estimates"))
    }else{
      pct <- readOGR(paste0("raw_data/vest/vest_2020/", tolower(state), "_2020"),
                     paste0(tolower(state), "_2020"))
    }
    pct <- spTransform(pct, CRS("+proj=longlat +datum=NAD83 +no_defs"))
    
    pct@data$GEOID <- paste0(state, c(1:nrow(pct@data)))
    pings  <- SpatialPoints(bgs@data[,c('INTPTLON10','INTPTLAT10')], proj4string = pct@proj4string)
    bgs$precinct <- over(pings, pct)$GEOID
    
    h20 <- bgs@data |> 
      mutate(state = state,
             year = 2020)
    
    saveRDS(bind_rows(h16, h20),
            paste0("temp/p2b_", state, ".rds"))
  #}
}

cl <- makeCluster(8)  
registerDoParallel(cl)

clusterExport(cl, list("b2p"))

states <- unique(filter(fips_codes, state_code <= "56")$state)

c(parLapply(cl, states,
            fun = b2p))


files <- list.files(path = "temp/", pattern = "^p2b*", full.names = T)

all_p2b <- rbindlist(lapply(files, readRDS)) |> 
  select(GEOID10, precinct, state, year)

saveRDS(all_p2b, "temp/precinct_block.rds")